Libraries and whatnot

library(tidyverse)
# library(broom)
library(patchwork)
library(scales)
library(dagitty)
library(ggdag)
library(latex2exp)  # Easily convert LaTeX into arcane plotmath expressions
library(ggtext)     # Use markdown in ggplot labels

# Create a cleaner serifed theme to use throughout
theme_do_calc <- function() {
  theme_dag(base_family = "Times New Roman") +
    theme(plot.title = element_text(size = rel(1.5)),
          plot.subtitle = element_markdown())
}

# Make all geom_dag_text() layers use these settings automatically
update_geom_defaults(ggdag:::GeomDagText, list(family = "Times New Roman", 
                                               fontface = "bold",
                                               color = "black"))

Examples

In all the examples that follow, we are trying to estimate a causal relationship between x and y. I think this is \(P(Y| \mbox{do}(X=x))\) NEED TO CHECK THIS

Chains (mediation)

This represents a sequential path from intervention to outcome via a confounder (x → z → y). Can bias causal link between x and y

Forks (common cause)

These occur when a confounding variable impacts both the intervention (x) and the outcome (Y). Can result in a spurious correlations (i.e. x and y will be correlated without existing causality).

Simple fork

This has no causal link between x and y, and a single confounder (z)

fork_dag <- dagify(x ~ z,
                   y ~ z,
                   coords = list(x = c(x = 1, y = 3, z = 2),
                                 y = c(x = 1, y = 1, z = 2)),
                   exposure = "x",
                   outcome = "y")

ggdag(fork_dag) +
  theme_dag()

Independence: x and y are independent (conditioned on z)

impliedConditionalIndependencies(fork_dag)
## x _||_ y | z

Adjustment set: shows what adjustment is necessary to obtain the causal link between input and outcome

ggdag_adjustment_set(fork_dag) +
  theme_dag()

Simulated data - Equations

  • z: \(N(10, 5)\)
  • x ~ z: \(\beta_1 = 5\)
  • y ~ z: \(\beta_2 = 2\)
n <- 1000
z <- rnorm(n, 10, 5)
x <- z * 5 + rnorm(n, 0, 2)
y <- z * 2 + rnorm(n, 0, 2)
df <- data.frame(x, y, z)
pairs(df)

Model of y ~ x. Ignoring the confounder results in spurious relationship between intervention and outcome

fit <- lm(y ~ x, df)
summary(fit)
## 
## Call:
## lm(formula = y ~ x, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6588 -1.4217  0.0332  1.4669  7.3531 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.036410   0.152095   0.239    0.811    
## x           0.399385   0.002756 144.920   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.172 on 998 degrees of freedom
## Multiple R-squared:  0.9546, Adjusted R-squared:  0.9546 
## F-statistic: 2.1e+04 on 1 and 998 DF,  p-value: < 2.2e-16

Model controlling for z. This removes the backdoor path between x and y and removes the apparent relationship

fit <- lm(y ~ x + z, df)
summary(fit)
## 
## Call:
## lm(formula = y ~ x + z, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.9982 -1.3117  0.0697  1.3173  5.9398 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.11968    0.14047  -0.852    0.394    
## x           -0.02330    0.03151  -0.739    0.460    
## z            2.12496    0.15790  13.458   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.999 on 997 degrees of freedom
## Multiple R-squared:  0.9616, Adjusted R-squared:  0.9615 
## F-statistic: 1.249e+04 on 2 and 997 DF,  p-value: < 2.2e-16

Simple fork 2

This is the same graph as before, but now we include a real relationship between x and y

fork_dag <- dagify(x ~ z,
                   y ~ x + z,
                   coords = list(x = c(x = 1, y = 3, z = 2),
                                 y = c(x = 1, y = 1, z = 2)),
                   exposure = "x",
                   outcome = "y")

ggdag(fork_dag) +
  theme_dag()

Independence: Now x and y have no conditional independence

impliedConditionalIndependencies(fork_dag)

Adjustment set:

ggdag_adjustment_set(fork_dag) +
  theme_dag()

Equations

  • z: \(N(10, 5)\)
  • x ~ z: \(\beta_1 = 5\)
  • y ~ x + z: \(\beta_2 = -4\), \(\beta_3 = 2\)
n <- 1000
z <- rnorm(n, 10, 5)
x <- z * 5 + rnorm(n, 0, 5)
y <- x * -4 + z * 2 + rnorm(n, 0, 5)
df <- data.frame(x, y, z)
pairs(df)

Model of y ~ x. Without controlling for the confounder, the coefficient relating x and y is biased (the relationship exists, but it’s biased)

fit <- lm(y ~ x, df)
summary(fit)
## 
## Call:
## lm(formula = y ~ x, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.6047  -3.8278   0.1458   3.8001  17.0101 
## 
## Coefficients:
##              Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)  0.360705   0.352867    1.022    0.307    
## x           -3.610218   0.006354 -568.171   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.316 on 998 degrees of freedom
## Multiple R-squared:  0.9969, Adjusted R-squared:  0.9969 
## F-statistic: 3.228e+05 on 1 and 998 DF,  p-value: < 2.2e-16

Model controlling for z. Now we get the right value for the coefficient (\(\sim4\))

fit <- lm(y ~ x + z, df)
summary(fit)
## 
## Call:
## lm(formula = y ~ x + z, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.6118  -3.2955  -0.0864   3.5797  15.0722 
## 
## Coefficients:
##             Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)  -0.3900     0.3343   -1.167    0.244    
## x            -3.9966     0.0318 -125.686   <2e-16 ***
## z             2.0040     0.1620   12.368   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.952 on 997 degrees of freedom
## Multiple R-squared:  0.9973, Adjusted R-squared:  0.9973 
## F-statistic: 1.861e+05 on 2 and 997 DF,  p-value: < 2.2e-16

Two variable fork

This is to illustrate Pearl’s point about uncessary confounders. This graph has a backdoor path with two confounding variables (z1 and z2). In order to test a causal link between x and y, we only need to control for one, as either will block the back door path

fork_dag <- dagify(x ~ z1,
                   z2 ~ z1,
                   y ~ z2,
                   coords = list(x = c(x = 1, y = 3, z1 = 1.5, z2 = 2.5),
                                 y = c(x = 1, y = 1, z1 = 2, z2 = 1.5)),
                   exposure = "x",
                   outcome = "y")

ggdag(fork_dag) +
  theme_dag()

Independence:

impliedConditionalIndependencies(fork_dag)
## x _||_ y | z2
## x _||_ y | z1
## x _||_ z2 | z1
## y _||_ z1 | z2

Adjustment set: Now we get two adjustment sets, as controlling for either z1 or z2 removes the path

ggdag_adjustment_set(fork_dag) +
  theme_dag()

Equations

  • z1: \(N(10, 5)\)
  • x ~ z1: \(\beta_1 = 5\)
  • y ~ z2: \(\beta_2 = 2\)
  • z2 ~ z1: \(beta_3 = -4\)
n <- 1000
z1 <- rnorm(n, 10, 5)
x <- z1 * 5 + rnorm(n, 0, 5)
z2 <- z1 * -4 + rnorm(n, 0, 5)
y <- z2 * 2 + rnorm(n, 0, 5)
df <- data.frame(x, y, z1, z2)
pairs(df)

Model of y ~ x (shows spurious relationship)

fit <- lm(y ~ x, df)
summary(fit)
## 
## Call:
## lm(formula = y ~ x, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -47.614  -9.764  -0.455   9.385  51.687 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.97585    0.97257   -3.06  0.00227 ** 
## x           -1.53535    0.01734  -88.53  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.8 on 998 degrees of freedom
## Multiple R-squared:  0.8871, Adjusted R-squared:  0.8869 
## F-statistic:  7838 on 1 and 998 DF,  p-value: < 2.2e-16

Control for z1

fit1 <- lm(y ~ x + z1, df)
summary(fit1)
## 
## Call:
## lm(formula = y ~ x + z1, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.088  -7.864   0.157   7.697  38.269 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.40923    0.83158   0.492    0.623    
## x           -0.04147    0.07402  -0.560    0.575    
## z1          -7.81078    0.37948 -20.583   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.56 on 997 degrees of freedom
## Multiple R-squared:  0.9207, Adjusted R-squared:  0.9206 
## F-statistic:  5790 on 2 and 997 DF,  p-value: < 2.2e-16

Control for z2

fit2 <- lm(y ~ x + z2, df)
summary(fit2)
## 
## Call:
## lm(formula = y ~ x + z2, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.8357  -3.7067   0.0417   3.5356  15.5855 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.04009    0.36113   0.111    0.912    
## x            0.01185    0.02049   0.579    0.563    
## z2           2.01444    0.02533  79.513   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.094 on 997 degrees of freedom
## Multiple R-squared:  0.9846, Adjusted R-squared:  0.9846 
## F-statistic: 3.19e+04 on 2 and 997 DF,  p-value: < 2.2e-16

Control for both (this works, but is more costly in regression terms)

fit3 <- lm(y ~ x + z1 + z2, df)
summary(fit3)
## 
## Call:
## lm(formula = y ~ x + z1 + z2, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.2111  -3.6316   0.0424   3.4860  15.0231 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.12908    0.36529  -0.353  0.72390    
## x           -0.05705    0.03251  -1.755  0.07956 .  
## z1           0.57546    0.21124   2.724  0.00656 ** 
## z2           2.06802    0.03201  64.604  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.078 on 996 degrees of freedom
## Multiple R-squared:  0.9847, Adjusted R-squared:  0.9847 
## F-statistic: 2.141e+04 on 3 and 996 DF,  p-value: < 2.2e-16

Colliders

This is the second type of path junction that can cause problems. Here, our two variables (x and y) are independent, but are both linked to z, the collider. Ths issue with this relationship is the bias that results when we control for the collider